module ground_control
  !
  ! Command Line Parser
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
#define green(i) (iand(mana,i)/=0)
  use suty
  implicit none
  private
  integer,public,protected::mana
  include "emily.f90"

  !
  public cmd_info
  public ppem_env
  public r4v3
  public me_first_mana
  public gc_veto
  public gc_veto_c
  public gc_veto_s
  public gc2mana
  public echo_v
  public echo_v3
  public echo_v4
  public etime
  public write_v_shm
  public write_v
  public fopen
  public get_einfo
  public set_rstat
  public mess_ioerr
  public proc_output
  public iostatus
  public nonzeronum

  !
  character*(*),public,parameter::ext_dcm=".dcm"
  character*(*),public,parameter::ext_raw=".img"

  !
  integer,parameter::GC_GREEN  =  0
  integer,parameter::GC_RED    = -1
  integer,parameter::GC_ORANGE =  1

  !
  integer,parameter,public::MAX_STR_LEN=512
  integer,parameter,public::MAX_LEN_FILENAME=512
  integer,parameter,public::MAXSTRLEN=MAX_STR_LEN
  integer,parameter,public::MAXPATHLEN=MAX_LEN_FILENAME

  !
  integer,parameter,public::FT_HIST=1
  integer,parameter,public::FT_LIST=2

  !
  interface echo_v
     module procedure echo_v3
     module procedure echo_v4
  end interface

  !
  type ppem_env
     real detector_separation
     integer num_crystal_x
     integer num_crystal_y
     integer num_crystal_x_in_detector
     integer num_crystal_in_gap
     integer num_detector_in_head
     integer num_doi_sampling
     real doi_sampling_step
     real voxel_height
     integer multiray_option
     integer sysmat_evaluator
     integer option
     integer vam ! voxel alignment mode
     real fa     ! filter amplitude
     integer maxcd ! maximum crystal difference
     real thj ! threshold for Cj
     real thi ! threshold for Ci
     character(MAX_LEN_FILENAME),pointer::fsysmat
     character(MAX_LEN_FILENAME),pointer::fhp,fhd,fht
     integer lo,up ! lower, upper slice with middle = 0
     real cf ! calibration factor
  end type ppem_env

  !
  type cmd_info
     character(MAX_STR_LEN*4) args ! <<<
     character(MAX_LEN_FILENAME),pointer::foutb ! base file name
     character(MAX_LEN_FILENAME),pointer::fnorm
     character(MAX_LEN_FILENAME),pointer::fcj,fci
     character(MAX_LEN_FILENAME),pointer::shm
     character(64) mc
     character(MAX_LEN_FILENAME)  ofmt
     character*4 oext
     integer,allocatable::argindxfile(:) ! holds filetype
     integer nfiles(2)
     integer narg
     integer num_iteration_max
     integer num_subsets
     integer ofmtid
     type(ppem_env) pem
  end type cmd_info

  !
  type r4v3
     real*4,pointer::a(:,:,:)
  end type r4v3

  !
  integer::ncall=0
  integer  rc
  character(MAX_STR_LEN) rep
  integer,parameter::ERR_INVALID_OPT   =  -1
  integer,parameter::ERR_INVALID_ARG   =  -2
  integer,parameter::ERR_MISSING_ARG   =  -3
  integer,parameter::ERR_INVALID_HIST  =  -4
  integer,parameter::ERR_INVALID_NORM  =  -5
  integer,parameter::ERR_INVALID_FOUT  =  -6
  integer,parameter::ERR_INVALID_INPUT =  -7
  integer,parameter::ERR_MISSING_SHM   =  -8
  integer,parameter::ERR_INVALID_SHM   =  -9
  integer,parameter::ERR_CPUID         = -10
  integer,parameter::ERR_MISSING_INPUT = -11
  integer,parameter::ERR_MISSING_OUTPUT= -12
  integer,parameter::ERR_INVALID_FORMAT= -13

  !
  integer,parameter::IOS_EEXIST = 17
  integer,parameter::IOS_EACCES = 13
  integer,parameter::IOS_NOENT  = 2

  !
  integer,parameter,public::EID_OPEN      =  1
  integer,parameter,public::EID_WRITE     =  2
  integer,parameter,public::EID_READ      =  3
  integer,parameter,public::EID_BUFOVER_1 = -1
  integer,parameter,public::EID_BUFOVER_2 = -2

  !
  integer,parameter,public::RSTAT_RUN    =   1
  integer,parameter,public::RSTAT_DONE   =   0
  integer,parameter,public::RSTAT_ERROR  = 999

  !
  integer mypid

  !
  type t_shmstat
     integer p
     integer max_size
     integer size
  end type t_shmstat
  integer::preconstat=0

  !
  integer,parameter::nshm=3
  integer,parameter,public::SHMID_ENV  = 3
  integer,parameter,public::SHMID_IN   = 1
  integer,parameter,public::SHMID_OUT  = 2
  type(t_shmstat),public::shmstats(nshm)

  !
  integer gcflag

  !
  integer,parameter,public::sp=selected_real_kind(6,37)
  integer,parameter,public::dp=selected_real_kind(15,307)

contains
  !
  subroutine print_info(to)
    integer,intent(in)::to
    write(*,1) PROGRAM//" "//VERSION
    if(to<=1) return
    write(*,1) COPYRIGHT//"\n"
    if(to<=2) return
    write(*,1) "Build at "//BUILD_TIME//" on "//HOST_NAME//"\n"
    if(to<=3) return
    write(*,1) NOTICE
1   format(a)
  end subroutine print_info
  !
  subroutine print_usage
    write(*,1) USAGE
1   format(a)
  end subroutine print_usage
  !
  character*12 function time_stamp()
    time_stamp=""
    call date_and_time(time=time_stamp)
    time_stamp="["//trim(time_stamp)//"]"
  end function time_stamp
  !
  subroutine echo(m,put_time)
    character*(*),intent(in)::m
    logical,intent(in),optional::put_time
    if(present(put_time).and.put_time) then
       write(*,20) time_stamp()
20     format(x,a,$)
    end if
    write(*,10) trim(m)
    call flush(6)
10  format(x,"(gc)",x,a)
  end subroutine echo
  !
  subroutine echo_red(msg,code)
    character*(*),intent(in),optional::msg
    integer,intent(in),optional::code
    integer ec
    character*(*),parameter::swan_song="*** mission aborted."
    if(present(msg)) call echo(trim(adjustl(msg)))
    call set_rstat(RSTAT_ERROR)
    write(0,"(x,a)") swan_song
    ec=-1
    if(present(code)) ec=code
    call exit(ec)
  end subroutine echo_red
  !
  subroutine proc_rangestr(str,ci,istat)
    character*(*),intent(in)::str
    type(cmd_info),intent(inout)::ci
    integer,intent(out)::istat
    integer k,lo,up
    istat=-1
    k=index(trim(str),":")
    if(k==0) then
       k=index(trim(str),",")
       if(k==0) return
    end if
    if(k==1.or.k==len_trim(str)) return
    read(str(1:k-1),*,iostat=istat) lo
    if(istat/=0) return
    read(str(k+1:),*,iostat=istat) up
    if(istat/=0) return
    if(lo>up) return
    istat=0
    ci%pem%lo=lo
    ci%pem%up=up
  end subroutine proc_rangestr 
  !
  subroutine init_shmstat
    shmstats(:)%p=0
    shmstats(:)%max_size=0
    shmstats(:)%size=0
  end subroutine init_shmstat
  !
  subroutine me_first_mana(ci)
    type(cmd_info),intent(inout),target::ci
    character(MAX_LEN_FILENAME),pointer::s
    real*8,pointer::x
    integer,pointer::k
    integer,pointer::in
    character(MAXSTRLEN) argv
    integer i,len,n,j
    integer filetype
    character(MAX_LEN_FILENAME),target::tmpstr
    external sproc
    pointer(psproc,sproc)
    !
    call init_shmstat
    call init_ci
    n=command_argument_count()
    ci%narg=n
    call get_command(ci%args)
    psproc=0
    mana=OPT_NOP
    mypid=getpid()
    rc=0
    rep=""
    filetype=FT_HIST ! histogram
    !
    nullify(s,x,k,in)
    allocate(ci%argindxfile(n))
    ci%argindxfile=0
    if(n==0) then
       call print_info(2)
       call print_usage
       call exit(0)
    end if
    main: do i=1,n
       call get_command_argument(i,argv,len)
       argv=trim(adjustl(argv))
       len=len_trim(argv)
       if(read_arg()) cycle
       if(rc/=0) exit
       if(argv(1:1)=="-".and.len>=2) then
          if(argv(2:2)/="-") then
             do j=2,len
                select case(argv(j:j))
                case("4")
                   mana = ior(mana, OPT_READ_32)
                case("a")
                   k => ci%pem%maxcd
                case("c")
                   mana = ior(mana, OPT_SHOW_STAT)
                case("d")
                   k=>ci%pem%num_detector_in_head
                case("g")
                   k=>ci%pem%num_crystal_in_gap
                case("f")
                   mana = ior(mana, OPT_ALLOW_OVERWRT)
                   if(green(OPT_RENAME)) call gc_veto_c(OPT_RENAME)
                case("h")
                   call print_usage
                   call exit(0)
                case("i")
                   in=>ci%argindxfile(i)
                   filetype=FT_HIST
                case("l")
                   in=>ci%argindxfile(i)
                   filetype=FT_LIST
                case("m")
                   k=>ci%pem%multiray_option
                case("n")
                   mana = ior(mana, OPT_NORMALIZE)
                   if(.not.associated(ci%fnorm)) allocate(ci%fnorm)
                   s=>ci%fnorm
                case("o")
                   allocate(ci%foutb)
                   s=>ci%foutb
                case("p")
                   x=>ci%pem%thj
                case("q")
                   call gc_veto_c(OPT_VERBOSE)
                case("s")
                   x=>ci%pem%doi_sampling_step
                case("t")
                   mana = ior(mana, OPT_DECAY_CORRECTION)
                case("u")
                   mana = ior(mana, OPT_CALIBRATION)
                case("v")
                   mana = ior(mana, OPT_VERBOSE)
                case("w")
                   mana = ior(mana, OPT_WRITE_IMAG)
                case("x")
                   k=>ci%pem%num_crystal_x
                case("y")
                   k=>ci%pem%num_crystal_y
                case("z")
                   k=>ci%pem%num_doi_sampling
                case("C")
                   k=>ci%pem%sysmat_evaluator
                case("D") 
                   x=>ci%pem%detector_separation
                case("E")
                   k=>ci%pem%vam
                case("F")
                   x=>ci%pem%fa
                   mana = ior(mana, OPT_USE_FILTER)
                case("G")
                   mana = ior(mana, OPT_NO_PAUSE)
                case("L")
                   mana = ior(mana, OPT_USE_LOGL)
                case("M")
                   mana = ior(mana, OPT_USE_SHM)
                   allocate(ci%shm)
                   s=>ci%shm                   
                case("N")
                   k=>ci%num_iteration_max
                case("O")
                   mana = ior(mana, OPT_USE_OSEM)
                case("P")
                   x=>ci%pem%thi
                case("R")
                   mana = ior(mana, OPT_RENAME)
                   if(green(OPT_ALLOW_OVERWRT)) call gc_veto_c(OPT_ALLOW_OVERWRT)
                case("S")
                   mana = ior(mana, OPT_SYSMAT_ONLY)
                case("T")
                   mana = ior(mana, OPT_USE_TRUE)
                case("U")
                   mana = ior(mana, OPT_UNIFORM_INIT)
                case("V")
                   call print_info(4)
                   call exit(0)
                case("W")
                   mana = ior(mana, OPT_WRITE_INTEGER)
                case("X")
                   k=>ci%pem%num_crystal_x_in_detector
                case("Z") 
                   x=>ci%pem%voxel_height
                case default
                   rc=ERR_INVALID_OPT
                   exit main
                end select
             end do
          else
             select case(argv(3:))
             case("Ci")
                if(.not.associated(ci%fci)) allocate(ci%fci)
                s=>ci%fci
                mana = ior(mana, OPT_LOAD_CI)
             case("Cj")
                if(.not.associated(ci%fcj)) allocate(ci%fcj)
                s=>ci%fcj
                mana = ior(mana, OPT_LOAD_CJ)
             case("cf")
                x => ci%pem%cf
             case("debug")
                mana = ior(mana, OPT_DEBUG)
             case("dump-Ci")
                if(.not.associated(ci%fci)) allocate(ci%fci)
                s=>ci%fci
                mana = ior(mana, OPT_DUMP_CI)
             case("dump-Ui")
                mana = ior(mana, OPT_DUMP_UI)
             case("dump-projection")
                mana = ior(mana, OPT_DUMP_PROJ) 
             case("flip-b1")
                mana = ior(mana, OPT_LMD_FLIPB1)
             case("flip-x1")
                mana = ior(mana, OPT_LMD_FLIPX1)
             case("flip-x2")
                mana = ior(mana, OPT_LMD_FLIPX2)
             case("flip-y1")
                mana = ior(mana, OPT_LMD_FLIPY1)               
             case("load-sysmat")
                if(.not.associated(ci%pem%fsysmat)) allocate(ci%pem%fsysmat)
                s=>ci%pem%fsysmat
                mana = ior(mana, OPT_LOAD_SYSMAT)
             case("ofmt")
                s => ci%ofmt
             case("range")
                s=>tmpstr
                psproc=loc(proc_rangestr)
             case("save-Cj")
                if(.not.associated(ci%fcj)) allocate(ci%fcj)
                s=>ci%fcj
                mana = ior(mana, OPT_SAVE_CJ)                
             case("save-Ci")
                if(.not.associated(ci%fci)) allocate(ci%fci)
                s=>ci%fci
                mana = ior(mana, OPT_SAVE_CI)
             case("save-sysmat")
                if(.not.associated(ci%pem%fsysmat)) allocate(ci%pem%fsysmat)
                s=>ci%pem%fsysmat
                mana = ior(mana, OPT_SAVE_SYSMAT)
             case("write-Cj")
                mana = ior(mana, OPT_WRITE_CJ)
             case("write-hist-true")
                allocate(ci%pem%fht)
                s=>ci%pem%fht
             case("write-hist-prompt")
                allocate(ci%pem%fhp)
                s=>ci%pem%fhp
             case("write-hist-delayed")
                allocate(ci%pem%fhd)
                s=>ci%pem%fhd
             case("write-update")
                mana = ior(mana, OPT_WRITE_UPDT)
             case default
                rc=ERR_INVALID_OPT
                exit main
             end select
          end if
       else
          ci%argindxfile(i)=filetype ! default is histogram
          ci%nfiles(filetype)=ci%nfiles(filetype)+1
       end if
    end do main

    if(green(OPT_VERBOSE)) write(*,*) trim(ci%args)

    !
    select case(rc)
    case(ERR_INVALID_OPT)
       rep="invalid option "//"'"//trim(argv)//"'"
    case(ERR_INVALID_ARG)
       rep="invalid argument "//"'"//trim(argv)//"'"
    case(GC_GREEN)
       if(associated(s).or.associated(x).or.associated(k)) then
          rep="missing argument for '"//trim(argv)//"'"
          rc=ERR_MISSING_ARG
       end if
    end select
    !
    call gc_veto(ci)
    !
  contains
    
    subroutine write_args
      integer ja
      character(MAXSTRLEN) a
      do ja=0,ci%narg
         call get_comamnd_argument(ja,a)
         write(*,1) trim(a)
1        format(x,a,$)
      end do
      write(*,*)
    end subroutine write_args
!!$    !
!!$    integer function get_narg()
!!$      integer pp
!!$      get_narg=0
!!$      do pp=1,len_trim(ci%args)
!!$         if(ci%args(pp:pp)==" ") get_narg=get_narg+1
!!$      end do
!!$    end function get_narg
!!$    !
!!$    subroutine get_arg(ii,aarg,llen)
!!$      integer,intent(in)::ii
!!$      character*(*),intent(out)::aarg
!!$      integer,intent(out)::llen
!!$      integer pp,nn,pp1,pp2
!!$      nn=1
!!$      pp1=1
!!$      pp2=0
!!$      do pp=1,len_trim(ci%args)
!!$         if(ci%args(pp:pp)==" ") then
!!$            nn=nn+1
!!$            if(ii==nn-1) then
!!$               pp2=pp-1
!!$               exit
!!$            else if(ii==nn) then
!!$               pp1=pp+1
!!$            end if            
!!$         end if
!!$      end do
!!$      if(pp2==0) then
!!$         llen=0
!!$         aarg=""
!!$      else
!!$         llen=pp2-pp1+1
!!$         aarg=ci%args(pp1:pp2)
!!$      end if
!!$    end subroutine get_arg
    !
    subroutine init_ci
      ci%nfiles=0
      ci%narg=0
      ci%args=""
      if(allocated(ci%argindxfile)) deallocate(ci%argindxfile)
      nullify(ci%fnorm)
      nullify(ci%foutb)
      nullify(ci%fcj)
      nullify(ci%fcj)
      nullify(ci%shm)
      ci%num_subsets=1
      ci%num_iteration_max=0
      ci%ofmtid=0
      ci%ofmt=valid_ofmt(1)
      ci%oext=""
      call init_pem(ci%pem)
    end subroutine init_ci
    !
    subroutine init_pem(pem)
      type(ppem_env),intent(out)::pem
      nullify(pem%fsysmat)
      nullify(pem%fhp)
      nullify(pem%fhd)
      nullify(pem%fht)
      pem%detector_separation=0
      pem%num_crystal_x=0
      pem%num_crystal_y=0
      pem%num_crystal_x_in_detector=0
      pem%num_crystal_in_gap=0
      pem%num_detector_in_head=0
      pem%num_doi_sampling=0
      pem%doi_sampling_step=0
      pem%voxel_height=1 ! <<<
      pem%multiray_option=0
      pem%sysmat_evaluator=0
      pem%option=0
      pem%fa=0
      pem%vam=0
      pem%maxcd=-1
      pem%thj=0
      pem%thi=0
      pem%lo=-huge(0)
      pem%up=huge(0)
      pem%cf=1
    end subroutine init_pem
    !
    logical function read_arg()
      integer iistat
      iistat=-1
      read_arg=.false.
      if(associated(s)) then
         s=trim(adjustl(argv))
         iistat=0
         if(psproc/=0) then
            call sproc(s,ci,iistat)
            psproc=0
         end if
         nullify(s)
      elseif(associated(x)) then
         read(argv,*,iostat=iistat) x
         nullify(x)
      elseif(associated(k)) then
         read(argv,*,iostat=iistat) k
         nullify(k)
      elseif(associated(in)) then
         ci%argindxfile(i)=filetype
         ci%nfiles(filetype)=ci%nfiles(filetype)+1
         nullify(in)
         iistat=0
      else
         return
      end if
      if(iistat==0) then 
         read_arg=.true.
      else
!         rc=ERR_INVALID_ARG THIS FAILS WITH -i or -l flags or file given
      end if
    end function read_arg
    !
  end subroutine me_first_mana
  !
  subroutine gc_veto(ci)
    use vuty
    type(cmd_info),intent(inout)::ci
    integer*4 buf(4)
    character(MAX_STR_LEN) errmsg
    character(MAXPATHLEN) f
    integer i,nf,ios
    !
    call gc2mana

    if(rc==0) then

       if(green(OPT_USE_OSEM)) ci%num_subsets=4
       if(ci%num_iteration_max==0.and..not.green(OPT_USE_LOGL)) then
          ci%num_iteration_max=1
          call echo("*** Warning number of iteration not set")
       end if
       if(green(OPT_RENAME).and.green(OPT_ALLOW_OVERWRT)) call gc_veto_c(OPT_ALLOW_OVERWRT)
    
       if(associated(ci%pem%fht).and..not.green(OPT_USE_TRUE)) call gc_veto_s(OPT_USE_TRUE)
          
       !
       if(green(OPT_LOAD_CJ).and.ci%pem%thi/=0.and.green(OPT_NORMALIZE)) &
            call gc_veto_c(OPT_LOAD_CJ) ! Cj must be calculated
       if(green(OPT_DEBUG)) then
          if(.not.green(OPT_VERBOSE)) call gc_veto_s(OPT_VERBOSE)
          if(.not.green(OPT_SHOW_STAT)) call gc_veto_s(OPT_SHOW_STAT)
       end if
       if(.not.green(OPT_NORMALIZE).and.green(OPT_LOAD_CI)) call gc_veto_c(OPT_LOAD_CI)
       if(green(OPT_LOAD_CJ)) call gc_veto_c(OPT_NORMALIZE)

       i=(ci%pem%num_crystal_x_in_detector+ci%pem%num_crystal_in_gap)*ci%pem%num_detector_in_head-ci%pem%num_crystal_in_gap
       if(ci%pem%num_crystal_x==0) then
          ci%pem%num_crystal_x=i
       else if(ci%nfiles(FT_LIST)>0) then
          if(ci%pem%num_crystal_x_in_detector==0.and.ci%pem%num_detector_in_head==0) then
             rep="*** Cannot determin num crystal x in detector."
             rc=ERR_MISSING_INPUT
          else if(ci%pem%num_crystal_x_in_detector==0) then
             if(mod(i-ci%pem%num_crystal_in_gap*(ci%pem%num_detector_in_head-1),&
                  ci%pem%num_detector_in_head)/=0) then
                rep="*** Cannot determin num crystal x in detector"
                rc=ERR_INVALID_INPUT
             else
                ci%pem%num_crystal_x_in_detector=(i-ci%pem%num_crystal_in_gap*(ci%pem%num_detector_in_head-1))&
                     /ci%pem%num_detector_in_head-ci%pem%num_crystal_in_gap                
             end if
          else if(ci%pem%num_detector_in_head==0) then
             if(mod(i+ci%pem%num_crystal_in_gap,ci%pem%num_crystal_x_in_detector+ci%pem%num_crystal_in_gap)/=0) then
                rep="*** Cannot determin num detector in head"
                rc=ERR_INVALID_INPUT
             else
                ci%pem%num_detector_in_head=(i+ci%pem%num_crystal_in_gap)&
                     /(ci%pem%num_crystal_x_in_detector+ci%pem%num_crystal_in_gap)
             end if
          else if(i/=ci%pem%num_crystal_x) then
             rep="*** Mismatch in num crystal x"
             rc=ERR_INVALID_INPUT             
          end if
       end if
       
       if(rc==0) then
          if(.not.check_range(ci%pem%detector_separation,0.0_dp,huge(1.0_dp),"detector separation")) then
             rep=errmsg
             rc=ERR_INVALID_INPUT
          else if(.not.check_range(real(ci%pem%num_crystal_x),0.0_dp,huge(1.0_dp),"num crystal x",4)) then
             rep=errmsg
             rc=ERR_INVALID_INPUT
          else if(.not.check_range(real(ci%pem%num_crystal_y),0.0_dp,huge(1.0_dp),"num crystal y",4)) then
             rep=errmsg
             rc=ERR_INVALID_INPUT
          else if(.not.check_range(real(ci%pem%thi),0.0_dp,1.0_dp,"threshold of Ci",1)) then
             rep=errmsg
             rc=ERR_INVALID_INPUT
          else if(.not.check_range(real(ci%pem%thj),0.0_dp,1.0_dp,"threshold of Cj",1)) then
             rep=errmsg
             rc=ERR_INVALID_INPUT
          else if(.not.green(OPT_SYSMAT_ONLY)) then
             if(green(OPT_USE_SHM)) then
                if(.not.check_shm(ci%shm)) then
                   rep="accessing shm: "//trim(errmsg)//", shmid="//trim(ci%shm)
                   rc=ERR_INVALID_SHM
                end if
             else
                nf=0
                do i=1,ci%narg
                   if(ci%argindxfile(i)>0) then
                      nf=nf+1
                      if(ci%argindxfile(i)==FT_HIST) then
                         call get_command_argument(i,f)
                         if(.not.check_fsize(f)) then
                            rep=trim(errmsg)//" histogram file "//trim(f)
                            rc=ERR_INVALID_HIST
                            exit
                         end if
                      end if
                   end if
                end do
                if(rc==0.and.nf==0) then
                   rc=ERR_MISSING_INPUT
                   rep=" missing histogram file"
                end if
             end if
          else
             rep=trim(rep)//" only proccess system matrix,"
          end if
       end if
    end if
    if(rc==0) then
       if(green(OPT_NORMALIZE).and..not.check_fsize(ci%fnorm)) then
          rc=ERR_INVALID_NORM
          rep=trim(errmsg)//" normalize file "//trim(ci%fnorm)
       end if
    end if
    if(rc==0) then
       if(associated(ci%foutb)) then
          do i=1,nofmt
             if(trim(ci%ofmt)==trim(valid_ofmt(i))) then
                ci%oext=ofmt_ext(i)
                ci%ofmtid=i
             end if
          end do
          if(ci%oext=="") then
             rc=ERR_INVALID_FORMAT
             rep="Invalid output format: "//trim(ci%ofmt)
          else
             f=ci%foutb
             if(proc_output(f,0,ios)==-1) then
                rc=ERR_INVALID_FOUT
                rep="failed to access output file "//trim(f)//" ["//trim(iostatus(ios))//"]"
             else
                rep=trim(rep)//" output "//trim(f)
                if(trim(f)/=trim(ci%foutb)) then
                   call echo("output file name altered: "//trim(ci%foutb)//" -> "//trim(f))
                   ci%foutb=f
                end if
             end if
          end if
       else if(green(OPT_SYSMAT_ONLY).or.green(OPT_USE_SHM)) then
          rep=trim(rep)//" no output file"
       else
          rc=ERR_MISSING_OUTPUT
          rep="missing output file"
       end if
    end if
    rep=adjustl(rep)
    call mana2gc
    
    !
    call gc2mana
    rep=""
    if(get_cpuid(buf,0)==GC_GREEN) then
       rc=get_cpuid(buf,int(z'80000002'))
       call int2str
       rc=get_cpuid(buf,int(z'80000003'))
       call int2str
       rc=get_cpuid(buf,int(z'80000004'))
       call int2str
       rep(len_trim(rep):)=""
    end if
    if(rc==GC_GREEN) then
       call mana2gc
    else
       rep="... couldnt get it started!"
       rc=ERR_CPUID
       call mana2gc
    end if
    !
    call gc2mana
    rep=""
    rc=get_cpuid(buf,1)
    rc=chex(buf(3),buf(4),ckdaz())
    call mana2gc
    !
    call gc2mana
    call regcode
    call mana2gc
    !
  contains    
    !
    subroutine regcode
      use kuty
      character*256 cksum
      integer*1 b(*)
      pointer(p,b)
      p=loc(ci%args)
      cksum=get_cksum(CKT_MD5,len_trim(ci%args),b(1:len_trim(ci%args)))
      ci%mc=trim(i2a(mana,hex=.true.))//":"//trim(cksum)
      rep="mission code "//trim(ci%mc)//" set"      
    end subroutine regcode
    !    
    logical function check_range(x,lo,up,name,type)
      real,intent(in)::x,lo,up
      character*(*),intent(in)::name
      integer,intent(in),optional::type ! 0: include bound, 1: exclude bound, 4: integer
      integer t
      character*32 fig
      check_range=.true.
      t=0
      if(present(type)) t=type
      
      if(iand(t,1)==0) then
         if(x<=lo) check_range=.false.
      else
         if(x<lo) check_range=.false.
      end if
      if(.not.check_range)then
         errmsg=trim(name)//" out of lower bounds"
         return
      end if

      if(iand(t,2)==0) then
         if(x>=up) check_range=.false.
      else
         if(x>up) check_range=.false.
      end if
      if(.not.check_range) then
         errmsg=trim(name)//" out of upper bounds"
         return
      end if
      rep=trim(rep)//" "//trim(name)
      if(iand(t,4)==0) then
         fig=f2a(x,-1)
      else
         fig=i2a(int(x))
      end if
      rep=trim(rep)//" "//trim(fig)//","
    end function check_range
    !
    logical function check_fsize(f)
      use kuty
      character*(*),intent(in)::f
      integer fs,hs
      hs=(ci%pem%num_crystal_x*ci%pem%num_crystal_y)**2
      if(green(OPT_READ_32)) then
         hs=hs*4
      else
         hs=hs*2
      end if
      check_fsize=.false.
      if(f=="") then
         errmsg="missing"
         return
      else
         fs=fsize(f)
         if(fs/=hs) then
            errmsg="size mismatch"
            return
         end if
         rep=trim(rep)//" "//trim(f)//" "//trim(i2a(fs))//","
      end if
      check_fsize=.true.
    end function check_fsize 
    !
    subroutine set_shmstat(p0,sid)
      use kuty
      integer,intent(in)::p0,sid
      integer p, qw
      pointer(pqw,qw)
      character*32 envstr
      character(MAXSTRLEN) str
      select case(sid)
      case(SHMID_OUT)
         envstr="RP3SHMOFFIMAG"
      case(SHMID_IN)
         envstr="RP3SHMOFFHIST"
      case(SHMID_ENV)
         envstr="RP3SHMOFFRECON"
      case default
         return
      end select
      shmstats(sid)%max_size=0
      shmstats(sid)%size=0
      shmstats(sid)%p=0
      call getenv(envstr,str)
      if(len_trim(str)>0) then
         p=a2i(str)
         if(p>=0) then
            p=p+p0
            shmstats(sid)%p=p
            pqw=p
            if(qw>0) then
               p=p+8
               shmstats(sid)%max_size=qw
               p=p+8
               shmstats(sid)%size=qw
            end if
         end if
      end if

    end subroutine set_shmstat
    !
    logical function check_shm(f)
      use kuty
      character*(*),intent(in)::f
      integer flg,id
      integer qw
      pointer(pqw,qw)
      integer hs,fs
      integer p0
      
      check_shm=.false.
      
      ! | int8 MAXSIZE | DATA |
      
      flg=int(o'0666',8)
      id=a2i(f)
      p0=shmat(id,0,0)
      if(p0<=0) then
         errmsg="shmat failed, error="//trim(i2a(geterrnoa()))
         return
      end if
      
      call set_shmstat(p0,SHMID_ENV)
      preconstat=shmstats(SHMID_ENV)%p
      call set_shmstat(p0,SHMID_OUT)
      call set_shmstat(p0,SHMID_IN)
      if(shmstats(SHMID_IN)%p>0) then
         fs=shmstats(SHMID_IN)%size
         if(fs>0) then
            hs=(ci%pem%num_crystal_x)**4
            if(green(OPT_READ_32)) then
               hs=hs*4
            else
               hs=hs*2
            end if
            if(hs/=fs) then
               errmsg="invalid size of histogram in shm sz="//trim(i2a(fs))
               return
            end if
         end if
      end if
      
      check_shm=.true.
    end function check_shm
    !
    subroutine int2str
      character*1 c
      integer*1 b
      integer i
      pointer(pc,c)
      pointer(pc,b)
      do i=0,15
         pc=loc(buf)+i
         if(b>=32) rep=trim(rep)//c
      end do
      rep=trim(rep)//"-"
    end subroutine int2str
    !
    integer function chex(ecx,edx,etc)
      integer*4,intent(in),target::ecx,edx,etc
      integer,target::edxbits(6),ecxbits(5),etcbits(1)
      character*8,target::edxlbl(6),ecxlbl(5),etclbl(1)
      integer,pointer::bits(:)
      integer*4,pointer::ex
      character*8,pointer::lbl(:)
      character*8 echo(2)
      integer ib
      integer veto
      integer k,i,j
      data edxlbl/"fpu","cmov","mmx","sse","sse2","htt"/
      data edxbits/0,15,23,25,26,28/
      data ecxlbl/"sse3","ssse3","sse41","sse42","popcnt"/
      data ecxbits/0,9,19,20,23/
      data etclbl/"daz"/
      data etcbits/0/
      data echo/"off","on"/
      veto=int(B'10000111111',8)
      chex=0
      bits=>edxbits
      lbl=>edxlbl
      ex=>edx
      ib=0
      do k=1,3
         do i=1,size(bits)
            rep=trim(rep)//" "//trim(lbl(i))
            if(btest(ex,bits(i))) then
               j=2
               veto=ibclr(veto,ib)
            else
               j=1
            end if
            rep=trim(rep)//" "//trim(echo(j))//", "
            ib=ib+1
         end do
         select case(k)
         case(1)
            bits=>ecxbits
            lbl=>ecxlbl
            ex=>ecx
         case(2)
            bits=>etcbits
            lbl=>etclbl
            ex=>etc
         end select
      end do
      chex=veto
      rep=adjustl(rep)
      rep(len_trim(rep):)=""
    end function chex
    !
    integer*4 function ckdaz()
      use vuty
      integer*4,allocatable::buf(:)
      integer istat
      ckdaz=0_4
      allocate(buf(128))
      buf=0_4
      call fxsave(loc(buf))
      if(btest(buf(8),6)) then
         ckdaz=1_4
         call setdaz(loc(istat))
      end if
      deallocate(buf)
    end function ckdaz
    !
  end subroutine gc_veto 
  !
  subroutine gc_veto_mod(i,set)
    integer,intent(in)::i
    logical,intent(in)::set
    character*64 str
    character*8 com
    integer k
    if(set) then
       k=ior(mana,i)
       com="SET"
    else
       k=iand(mana,not(i))
       com="CLEAR"
    end if
    if(k/=mana) then
       str="*** WARNING, GC-VETO "//trim(com)//" CODE, "//trim(i2a(i,hex=.true.))
       call echo(str)
       mana = k
    end if
  end subroutine gc_veto_mod
  !
  subroutine gc_veto_c(i)
    integer,intent(in)::i
    call gc_veto_mod(i,.false.)
  end subroutine gc_veto_c
  !
  subroutine gc_veto_s(i)
    integer,intent(in)::i
    call gc_veto_mod(i,.true.)
  end subroutine gc_veto_s
  !
  subroutine gc2mana(anc)
    character*(*),intent(in),optional::anc
    integer,parameter::npre=9
    character(MAX_STR_LEN) m
    character*128 hdr
    integer trem
    logical  pause
    character*1 ans
    character(lenbc) base
    pointer(pbase,base)
    character(lensc) ship
    pointer(pship,ship)
    character(lenpc) pref
    integer ls
    integer,parameter::lbye=27
    integer*1::bye(lbye)=[97_1,110_1,100_1,32_1,109_1,97_1,121_1,32_1,103_1,111_1,&
         100_1,115_1,32_1,108_1,111_1,118_1,101_1,32_1,98_1,101_1,&
         32_1,119_1,105_1,116_1,104_1,32_1,117_1]
    character(lbye) b
    pointer(ppref,pref)
    pointer(p,b)
    ncall=ncall+1
    if(ncall<=NMSGMAX) then
       if(green(OPT_DEBUG)) then
          pbase=loc(basecode)
          if(ncall==1) then
             pship=loc(shipcode)
             ls=lensc
          else
             pship=loc(captaincode)
             ls=lencc
          end if
          hdr=base//" to "//ship(1:ls)
          if(ncall==1) then
             ppref=loc(prefcode)
             hdr=pref//" "//trim(hdr)
          end if
          call echo(trim(hdr)//", "//trim(themsg(ncall,MSG_GC)),ncall==1)
       end if
    else
       trem=(npre+NMSGMAX+1)-ncall
       if(green(OPT_VERBOSE)) then
          pause=.not.green(OPT_NO_PAUSE).and.green(OPT_DEBUG)
          if(trem==0) then
             m="lift off."
          else
             if(trem<0) then
                m="Iteration "//trim(i2a(abs(trem)+1))//" begin..."
             else
                m=i2a(abs(trem))
                if(present(anc)) m=trim(m)//"... "//trim(anc)
                if(trem==1) then
                   p=loc(bye)
                   m=trim(m)//", "//b
                end if
                m=trim(m)//"."
             end if            
          end if
          if(pause) write(*,*)
          call echo(m, .true.)
          if(pause) then
             ans=trim(getans("ready? [Quit/Go/or Just Hit Enter] > "))
          else
             ans=""
          end if
          select case(ans)
          case("q","Q")
             call echo_red
          case("g","G")
             call gc_veto_s(OPT_NO_PAUSE)
          end select
       end if
    end if
  end subroutine gc2mana
  !
  subroutine mana2gc
    character*128 hdr
    character(MAX_STR_LEN) m
    character(lencc) cap
    character(lenbc) base
    character(lenpc) pref
    pointer(pcap,cap)
    pointer(pbase,base)
    pointer(ppref,pref)
    if(ncall>NMSGMAX) return
    m=themsg(ncall,MSG_CAP)
    pcap=loc(captaincode)
    pbase=loc(basecode)
    hdr=cap//" to "//base
    if(rep/="") m=trim(m)//" "//trim(rep)
    if(rc==GC_GREEN) then
       m=trim(m)//", green."
    else
       m=trim(m)//"."
    end if
    m=trim(hdr)//", "//trim(m)
    if(ncall==1) then
       ppref=loc(prefcode)
       m=pref//" "//trim(m)
    end if
    if(green(OPT_DEBUG).or.rc<=GC_RED) &
         call echo(m,ncall==1)
    if(rc<=GC_RED) call echo_red()
  end subroutine mana2gc
  !
  character*128 function themsg(i,gcap)
    integer,intent(in)::i,gcap
    integer j,k
    integer*1 b
    character*1 c
    pointer(p,b)
    pointer(p,c)
    themsg=""
    j=i
    p=0
    select case(gcap)
    case(MSG_GC)
       p=loc(gcmsg)
    case(MSG_CAP)
       p=loc(capmsg)
    end select
    p=p-1
    do while(j>=2)
       p=p+1
       if(b==-1) return
       if(b==0) j=j-1
    end do
    k=0
    do 
       p=p+1
       if(b==0) return
       k=k+1
       themsg(k:k)=c
    end do
  end function themsg
  !
  real*4 function nonzeromin(pv,nv)
    integer,intent(in),value::pv,nv
    integer i
    real*4 v
    pointer(p,v)
    p=pv-4
    gcflag=0
    nonzeromin=real(huge(0.0),4)
    do i=1,nv
       p=p+4
       if(v<nonzeromin) then
          if(v/=0) then
             nonzeromin=v
          else
             gcflag=1 ! zero found
          end if
       end if
    end do
  end function nonzeromin
  !
  integer function nonzeronum(pv,nv)
    integer,intent(in),value::pv,nv
    integer i
    real*4 v
    pointer(p,v)
    p=pv-4
    nonzeronum=0
    do i=1,nv
       p=p+4
       if(v/=0) nonzeronum=nonzeronum+1
    end do
  end function nonzeronum
  !
  subroutine echo_v3(v,t)
    real*4,intent(in)::v(:,:,:)
    character*(*),intent(in)::t
    if(.not.green(OPT_SHOW_STAT)) return
    call echo_vstat(nonzeromin(loc(v),sizeof(v)/4),maxval(v),sum(v),&
         nonzeronum(loc(v),sizeof(v)/4),t)
  end subroutine echo_v3
  subroutine echo_v4(v,t)
    real*4,intent(in)::v(:,:,:,:)
    character*(*),intent(in)::t
    if(.not.green(OPT_SHOW_STAT)) return
    call echo_vstat(nonzeromin(loc(v),sizeof(v)/4),maxval(v),sum(v),&
         nonzeronum(loc(v),sizeof(v)/4),t)
  end subroutine echo_v4
  subroutine echo_vstat(vmin,vmax,vsum,num,t)      
    real*4,intent(in)::vmin,vmax,vsum
    integer,intent(in)::num
    character*(*),intent(in)::t
    character*32 note
    if(gcflag==1.and.vmin>0) then
       ! zero found and all positive
       note="[mingt0/max/avegt0/sum]:"
    else
       !
       note="[min/max/avegt0/sum]:"
    end if
    write(*,10) trim(t)//" "//trim(note),vmin,vmax,vsum/num,vsum
    call flush(6)
10  format(x,a,4(x,es12.4))
  end subroutine echo_vstat
  !
  real function etime(t0,t1)
    integer,intent(in)::t0(:),t1(:)
    integer t(8)
    t=t1-t0
    if(t(8)<0) then
       t(8)=t(8)+1000
       t(7)=t(7)-1
    end if
    if(t(7)<0) then
       t(7)=t(7)+60
       t(6)=t(6)-1
    end if
    if(t(6)<0) then
       t(6)=t(6)+60
       t(5)=t(5)-1       
    end if
    if(t(5)<0) then
       t(5)=t(5)+24
       t(3)=t(3)-1
    end if
    etime=(t(7)+real(t(8))*1.e-3)+60*(t(6)+60*(t(5)+24*t(3)))
  end function etime
  !
  subroutine set_fn(f,e,i1,i2)
    character*(*),intent(inout)::f
    character*(*),intent(in)::e
    integer,intent(in)::i1,i2
    integer i
    if(i1>=0) then
       f=trim(f)//"-"//trim(i2a(i1))
    end if
    if(i2>0) then
       f=trim(f)//"_"//trim(i2a(i2))
    end if
    i=len_trim(f)
    if(f(i-3:i)/=e(1:4)) f=trim(f)//trim(e)
  end subroutine set_fn
  !
  subroutine write_v_shm(pvec,ifo)
    use duty
    use suty
    integer,intent(in)::pvec
    type(image_info),intent(in)::ifo
    integer*8 fs
    integer p
    integer*4 dw
    real*4 st
    p=shmstats(SHMID_OUT)%p+8

    ! | int8 max_size | int8 size | int4 width | int4 height | int4 number_of_slices |
    ! | real4 dx | real4 dy | real4 slice_width | real4 0 | real4 data |

    fs=ifo%height*ifo%width*ifo%number_of_slices*4
    p=mcpf(p,loc(fs),8)

    dw=ifo%width
    p=mcpf(p,loc(dw),4)
    dw=ifo%height
    p=mcpf(p,loc(dw),4)
    dw=int(ifo%number_of_slices,4)
    p=mcpf(p,loc(dw),4)
    dw=0_4
    p=mcpf(p,loc(dw),4)

    st=real(ifo%dx,4)
    p=mcpf(p,loc(st),4)
    st=real(ifo%dy,4)
    p=mcpf(p,loc(st),4)
    st=real(ifo%slice_width,4)
    p=mcpf(p,loc(st),4)
    st=0_4
    p=mcpf(p,loc(st),4)
    
    p=mcpf(p,pvec,fs)
  end subroutine write_v_shm
  !
  subroutine write_v(f,i1,i2,pvec,nvec,istat,fmtid,ifo)
    use duty
    use suty
    use auty
    character*(*),intent(in)::f
    integer,intent(in)::i1,i2
    integer,intent(in)::pvec,nvec
    integer,intent(out)::istat
    integer,intent(in)::fmtid
    type(image_info),intent(inout),optional::ifo
    real,parameter::int_max=16383
    real,parameter::int_ave=8192
    character(MAXPATHLEN) fout,fhdr,str
    integer*2,allocatable::iv(:)
    character*4 e
    real*4 v(*)
    real::norm=0,maxv
    logical dicom
    integer iu
    pointer(pv,v)
    type(ana_hdr) hdr
    integer i
    pv=pvec
    istat=0
    dicom=(fmtid==OFMT_DCM).and.present(ifo)
    fout=f
    e=ofmt_ext(fmtid)
    if(index(fout,".")/=0) e="" ! only one ext allowed
    call set_fn(fout,e,i1,i2)
    iu=proc_output(fout,10,istat)
    if(iu>0) then
       if(green(OPT_WRITE_INTEGER).or.dicom) then
          allocate(iv(nvec))
          maxv=maxval(v(1:nvec))
          if(maxv>0) norm=int_max/maxv
          iv(1:nvec)=int(v(1:nvec)*norm+0.5_dp,2)
          if(dicom) then
             ifo%intcpt=0
             ifo%bias=1.0_dp/norm
             istat=write_dcm(iu,ifo,loc(iv),nvec)
             call echo("Dicom integerized by: "//trim(f2a(norm,-1)))
          else
             write(iu,iostat=istat) iv(1:nvec)             
          end if
          deallocate(iv)
       else 
          write(iu,iostat=istat) v(1:nvec)
       end if
       close(iu)
       if(fmtid==OFMT_ANA.and.present(ifo)) then
          call msets(loc(hdr),0,sizeof(hdr))
          hdr%hk%data_type="dsr"
          hdr%hk%db_name=trim(fout)
          hdr%hk%sizeof_hdr=348_4
          hdr%hk%extents=16384_4
          hdr%hk%regular="r"
          hdr%id%dim(1)=4_2
          hdr%id%dim(2)=int(ifo%width,2)
          hdr%id%dim(3)=int(ifo%height,2)
          hdr%id%dim(4)=int(ifo%number_of_slices,2)
          hdr%id%dim(5)=1_2
          hdr%id%pixdim(1)=4_4
          hdr%id%pixdim(2)=real(ifo%dx,4)
          hdr%id%pixdim(3)=real(ifo%dy,4)
          hdr%id%pixdim(4)=real(ifo%slice_width,4)
          hdr%id%vox_units="mm"
          hdr%id%cal_units="count"
          hdr%id%datatype=int(DT_FLOAT,2)
          hdr%id%bitpix=32_2
          hdr%id%cal_max=1_4
          hdr%id%glmax=1_4
          call getenv("RP3STUDYDESC",hdr%dh%descrip)
          call getenv("RP3DERICDESC",hdr%dh%generated)
          hdr%dh%originator="emily"
          call getenv("RP3SID",hdr%dh%scannum)
          call getenv("RP3PID",hdr%dh%patient_id)
          call getenv("RP3SYEAR",str)
          hdr%dh%exp_date=trim(str)
          call getenv("RP3SMONTH",str)
          hdr%dh%exp_date=trim(hdr%dh%exp_date)//trim(str)
          call getenv("RP3SDAY",str)
          hdr%dh%exp_date=trim(hdr%dh%exp_date)//trim(str)
          call getenv("RP3SHOUR",str)
          hdr%dh%exp_time=trim(str)
          call getenv("RP3SMINUTE",str)
          hdr%dh%exp_time=trim(hdr%dh%exp_time)//trim(str)
          call getenv("RP3SSECOND",str)
          hdr%dh%exp_time=trim(hdr%dh%exp_time)//trim(str)       
          i=len_trim(fout)
          if(i>4.and.fout(i-3:i)==".img") then
             fhdr=fout(1:i-4)//".hdr"
          else
             fhdr=trim(fout)//".hdr"
          end if
          open(unit=2,file=fhdr,access="stream",iostat=istat,status="replace")
          if(istat/=0) call echo_red( "*** Error during open: " //trim(fhdr))
          write(2,iostat=istat) hdr
          if(istat/=0) call echo_red( "*** Error during write: " //trim(fhdr))
          close(2)
       end if
    end if
    call echo("Writing: "//trim(fout)//" ["//trim(iostatus(istat))//"]")
    if(istat/=0) call echo_red("write_v failed")
  end subroutine write_v
  !
  character*32 function iostatus(istat)
    integer,intent(in)::istat
    select case(istat)
    case(0)
       iostatus="ok"
    case default
       write(iostatus,*) istat
       iostatus="error code = "//trim(adjustl(iostatus))
    end select
  end function iostatus
  !
  subroutine fopen(f,i)
    character*(*),intent(in)::f
    integer,intent(in)::i
    integer istat
     open(unit=i,file=f,status="replace",iostat=istat)
     if(istat/=0) then
        call echo_red("*** cannot open file: "//trim(f))
     end if
   end subroutine fopen
   !
   subroutine get_einfo(ifo)
     use duty
     type(image_info),intent(inout)::ifo
     character*512 str
     call getenv("RP3PID",ifo%pinf%id)
     call getenv("RP3PNAME",ifo%pinf%name)
     call getenv("RP3PDOB",ifo%pinf%dob)
     call getenv("RP3PSEX",ifo%pinf%sex)
     call getenv("RP3PAGE",str)
     ifo%pinf%age=a2i(str)
     call getenv("RP3PWEIGHT",str)
     ifo%pinf%weight=a2i(str)
     call getenv("RP3PHEIGHT",str)
     ifo%pinf%height=a2i(str)
     call getenv("RP3IMGORIPAT",ifo%pinf%orientation)
     !
     call getenv("RP3SID",ifo%sinf%id)
     ifo%sinf%input=trim(ifo%sinf%id) ! <<<<<<<<<<<
     ifo%sinf%ofile=trim(ifo%sinf%input)//".dcm" ! <<<<<
     call getenv("RP3SYEAR",str)
     ifo%sinf%year=a2i(str)
     call getenv("RP3SMONTH",str)
     ifo%sinf%month=a2i(str)
     call getenv("RP3SDAY",str)
     ifo%sinf%day=a2i(str)
     call getenv("RP3SHOUR",str)
     ifo%sinf%hour=a2i(str)
     call getenv("RP3SMINUTE",str)
     ifo%sinf%minute=a2i(str)
     call getenv("RP3SSECOND",str)
     ifo%sinf%second=a2i(str)

     call getenv("RP3SPART",ifo%sinf%part)
     call getenv("RP3STUDYDESC",ifo%sinf%description)
!!$     if(len_trim(ifo%sinf%part)>0.and.len_trim(ifo%sinf%description)>0)&
!!$          ifo%sinf%description=trim(ifo%sinf%description)//" "//trim(ifo%sinf%part)
     call getenv("RP3DERIVDESC",ifo%einf%derivDesc)
     call getenv("RP3SINT",ifo%sinf%injectionTime)
     call getenv("RP3SDOSE",str)
     ifo%sinf%dose=a2f(str)
     call getenv("RP3SSCAND",str)
     ifo%sinf%duration=a2i(str)
     call getenv("RP3RPHAR",ifo%sinf%rphar)
     call getenv("RP3HALFLIFE",str)
     ifo%sinf%halflife=a2f(str)*60 ! in sec
     call getenv("RP3PFRAC",str)
     ifo%sinf%positronFrac=a2f(str) ! in sec
     !
     call getenv("RP3INSTITUTION",ifo%einf%institution)
     call getenv("RP3OPERATOR",ifo%einf%operator)
   end subroutine get_einfo

  !
   subroutine set_rstat(code)
     integer,intent(in)::code
     integer q
     pointer(p,q)
     if(preconstat==0) return
     p=preconstat
     q=code
   end subroutine set_rstat

   !
   character*128 function mess_ioerr(eid,f,i)
     integer,intent(in)::eid
     character*(*),intent(in)::f
     integer,intent(in)::i
     select case(eid)
     case(EID_OPEN)
        mess_ioerr="OPEN"
     case(EID_WRITE)
        mess_ioerr="WRITE"
     case(EID_READ)
        mess_ioerr="READ"
     end select
     mess_ioerr="*** Error during "//trim(mess_ioerr)//": "&
          //trim(f)//", error="//trim(i2a(i))
   end function mess_ioerr

   integer function test_output(f)
     character*(*),intent(in)::f
     integer istat
     open(unit=100,file=f,status="new",access="stream",iostat=istat)
     close(100)
     if(istat==0) call unlink(f)
     test_output=istat
   end function test_output

   integer function proc_output(f,iu,istat,ext)
     ! -1: failure
     !  0: success and not opened
     ! 0<: success and opened connected to the unit
     character*(*),intent(inout)::f
     integer,intent(in)::iu
     integer,intent(out),optional::istat
     character*(*),intent(in),optional::ext
     integer istat_
     logical batch
     proc_output=-1
     if(present(istat)) istat=0
     if(present(ext)) call putext
     do
        istat_= test_output(f)
        if(istat_==IOS_EEXIST) then
           call echo("*** File already exists: "//trim(f))
           if(.not.green(OPT_ALLOW_OVERWRT)) then
              batch=.not.green(OPT_VERBOSE)
              if(.not.batch) then
                 if(.not.green(OPT_RENAME)) then
                    select case(getans("Cancel/Retry/Overwrite/New/Auto? => "))
                    case("C","c")
                       call echo_red("proc_output: canceled by the user")
                    case("O","o")
                       call gc_veto_s(OPT_ALLOW_OVERWRT) ! <<< for all files!
                       istat_=0
                       exit
                    case("R","r")
                       cycle
                    case("N","n")
                       f=getans("File name => ")
                       cycle
                    case("A","a")                 
                       batch=.true.
                    case default
                       cycle
                    end select
                 end if
              end if
              call suggest_another(f)
           else
              istat_=0
              exit
           end if
        else
           exit
        end if
     end do

     if(iu>0.and.istat_==0) then
        open(unit=iu,file=f,status="replace",access="stream",iostat=istat_)
     end if
     
     if(istat_/=0) then
        proc_output=-1
     else
        if(iu>0) then
           proc_output=iu
        else
           proc_output=0
        end if
     end if
     if(present(istat)) istat=istat_

   contains
     
     subroutine putext
       integer lf,le
       lf=len_trim(f)
       le=len_trim(ext)
       if(lf<=le.or.f(lf-le+1:lf)/=ext(1:le)) then
          f=trim(f)//trim(ext)
       end if
     end subroutine putext

     subroutine suggest_another(fn)
       character*(*),intent(inout)::fn
       integer i,k,j
       character*(MAXPATHLEN) cext
       integer*1 b
       integer indx,lindx,ln,lf
       character*32 cindx
       pointer(p,b)
       k=0
       lf=len_trim(fn)
       k=lf+1
       cext=""
       do i=k,2,-1
          if(fn(i:i)==".") then
             k=i
             cext=fn(i:lf)
             exit
          end if
       end do
       ln=0
       j=k
       do i=k-1,2,-1
          p=loc(fn(i:i))
          if(b>=z"30".and.b<=z"39") then
             ln=ln+1
          else 
             if(b/=z"3A") then
                ln=0
             else
                j=i
             end if
             exit            
          end if                   
       end do
       if(ln>0) then
          indx=a2i(fn(j+1:j+ln))+1
       else
          indx=1
       end if
       cindx=":"//trim(i2a(indx))
       lindx=len_trim(cindx)
       fn(j:j+lindx-1)=trim(cindx)
       fn(j+lindx:)=trim(cext)
     end subroutine suggest_another

   end function proc_output
   
   character*(MAXSTRLEN) function getans(m)
     character*(*),intent(in)::m
     write(*,10) m
     read(*,20) getans
     10 format(a,$)
     20 format(a)
   end function getans
   
 end module ground_control
